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Abstract — Generalized  Fourier  series  with  orthogonal 
polynomial  bases  have  useful  applications  in  several  fields, 
including  differential  equations,  pattern  recognition,  and 
image  and  signal  processing.  However,  computing  the 
generalized  Fourier  series  can  be  a  challenging  problem, 
even  for  relatively  well  behaved  functions.  In  this  paper,  a 
method  for  approximating  a  sparse  collection  of  Fourier- 
like  coefficients  is  presented  that  uses  a  collocation  tech¬ 
nique  combined  with  an  optimization  problem  inspired  by 
recent  results  in  compressed  sensing  research.  The  dis¬ 
cussion  includes  approximation  error  rates  and  numerical 
examples  to  illustrate  the  effectiveness  of  the  method.  One 
example  displays  the  accuracy  of  the  generalized  Fourier 
series  approximation  for  several  test  functions,  while  the 
other  is  an  application  of  the  generalized  Fourier  series 
approximation  to  rotation-invariant  pattern  recognition  in 
images. 

I.  Introduction 

This  paper  discusses  an  efficient  method  to  approx¬ 
imate  the  sparse  generalized  Fourier  series  of  a  given 
function  in  terms  of  orthogonal  polynomials  in  several 
dimensions.  That  is,  given  a  function  /  :  — >  M 

satisfying  certain  properties,  we  seek  to  compute  the 
Fourier-like  coefficients  {cn  :  n  G  W}  such  that 

/  ~  y2  Cnvrn,  (1) 

n£W 

where  {7rn}  is  a  collection  of  orthogonal  polynomials 
and  W  is  is  a  sparse  collection  of  multi-indices. 

To  compute  the  sparse  collection  of  Fourier-like  co¬ 
efficients,  we  propose  using  a  collocation  method  in 
an  optimization  scheme  motivated  by  recent  results  in 
the  field  of  compressed  sensing.  This  method  seeks  a 
sparse  collection  of  coefficients  such  that  the  finite  sum 
appealing  in  (1)  nearly  interpolates  /  at  carefully  chosen 
nodes. 

The  notion  of  sparsity  has  been  increasingly  em¬ 
phasized  in  recent  literature.  As  the  amount  of  useful 
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data  collected  continues  to  increase,  the  extraction  of 
meaningful  data  is  often  hindered  by  the  so-called  ‘curse 
of  dimensionality’.  Too  alleviate  its  effects,  methods 
exploiting  any  natural  sparse  or  hierarchical  structures 
in  data  have  been  developed.  Research  in  sparse  grids, 
particularly  when  coupled  with  a  hierarchical  scheme  us¬ 
ing  Smolyak’s  algorithm,  laid  the  foundation  for  accurate 
computation  of  high  dimensional  problems  while  sig¬ 
nificantly  reducing  the  computational  complexity  from 
0(Nd)  to  0(N  log^1  N),  where  d  is  the  underlying 
dimensionality  of  the  problem  and  N  is  the  number  of 
grid  points  in  each  coordinate  direction  [3],  [19]. 

For  lineal-  regression  problems,  the  importance  of 
sparsity  has  emerged  in  recent  advances  in  compressed 
sensing.  Given  a  sensing  matrix  X  e  Wnxp  with  m  <C  p 
and  a  collection  of  under  sampled  linear  measurements 
y  =  Xc,  using  compressed  sensing  techniques  one  may 
approximate  the  sparse  vector  c  via  the  optimization 
problem 

c  =  argmin  ||c||i  satisfying  Xc  =  y. 

Under  certain  sparsity  and  coherence  conditions,  the 
recovery  will  be  exact  [2],  [4],  [5],  [7]. 

On  the  other  hand,  sparse  generalized  Fourier  series 
approximations  with  orthogonal  polynomial  bases  typi¬ 
cally  use  spectral  or  pseudospectral  approximations  on 
sparse  grids  [6],  [9],  [11],  [15],  [16].  These  methods 
determine  the  sparse  collection  of  multi-indices  a  priori, 
which  may  lead  to  unintentional  deletion  of  significant 
terms  in  the  approximation.  To  compute  the  coeffi¬ 
cients  in  the  aforementioned  works,  detailed  hierarchical 
quadrature  methods  have  been  developed  to  perform 
the  numerical  integration  required  for  the  desired  accu¬ 
racy  while  keeping  the  computational  complexity  low. 
Although  effective,  these  methods  can  be  difficult  to 
implement. 

In  comparison,  the  method  proposed  in  this  paper  has 
a  straightforward  implementation.  In  addition,  a  major 
strength  of  the  proposed  method  is  that  the  sparse  index 
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set  does  not  need  to  be  selected  in  advance.  Instead,  the 
optimization  scheme  will  select  the  best  sparse  selection 
of  multi-indices  to  approximate  the  function. 

The  following  notation  will  be  used  throughout  this 
work.  Vectors  will  be  denoted  by  boldface  letters.  Let 
M  denote  the  real  numbers,  N  the  natural  numbers  and 
No  =  NU{0}.  Let  Rd  denote  the  usual  d-dimensional  Eu¬ 
clidean  space  with  obvious  analogues  Nd  and  Nq.  Let  rij 
denote  the  jth  entry  of  the  vector  n.  Let  C(Fl)  equal  the 
space  of  all  continuous  real-valued  functions  defined  on 
a  domain  fh  The  space  of  unweighted  square-Lebesgue 
integrable  functions  over  is  denoted  by  (W1)  with 
the  associated  weighted  inner  product  (•,  -)u  defined 
by (f,g)u  =  Juifg  du;  for  all  f.g  £  Ll(Rd).  The  u- 
weighted  inner  product  naturally  induces  the  weighted 
norm  ||  •  ||w. 

The  rest  of  this  paper  is  organized  as  follows.  In 
Section  2,  the  proposed  method  for  approximating  the 
sparse  generalized  Fourier  series  of  a  function  is  pre¬ 
sented,  along  with  an  analysis  of  the  approximation 
error.  The  advantages  of  the  proposed  method  are  stated 
and  compared  to  the  difficulties  exhibited  by  existing 
methods.  Section  3  includes  numerical  experiments  in 
which  the  proposed  method  is  implemented.  The  first 
set  of  experiments  approximates  the  sparse  generalized 
Fourier  series  of  several  bivariate  test  functions.  The  sec¬ 
ond  set  of  experiments  applies  the  sparse  approximated 
Fourier-like  coefficients  to  compute  the  rotation-invariant 
Gaussian-Hermite  moments,  and  uses  these  moments 
to  perform  classification  of  unknown  rotated  images. 
The  paper  concludes  with  a  summary  and  directions  of 
possible  future  work. 


II.  Sparse  Generalized  Fourier  Series  with 
Orthogonal  Polynomial  Bases 

The  proposed  method  approximates  a  given  function 
by  a  sparse  generalized  Fourier  series  with  an  orthog¬ 
onal  polynomial  basis.  After  selecting  a  set  of  multi¬ 
indices  and  evaluating  the  orthogonal  polynomials  at  pre¬ 
determined  nodes,  the  sparse  collection  of  coefficients  is 
computed  using  a  collocation  method  with  a  convex  opti¬ 
mization  problem.  This  section  explores  the  multi-index 
and  node  selection,  thoroughly  explains  the  collocation 
model  used  to  approximate  the  Fourier-like  coefficients 
and  discusses  the  approximation  errors.  We  begin  with 
preliminary  facts  about  orthogonal  polynomials  and  gen¬ 
eralized  Fourier  series. 


A.  Full  and  sparse  grid  orthogonal  polynomial  repre¬ 
sentations 

Let  C  Rd.  Let  uj  :  fl  — >  M  be  a  continuous  weight 
function.  That  is,  suppose  u  >  0  on  Fl  and  uj  satisfies 


Associated  with  each  weight  function  is  a  collection  of 
orthogonal  polynomials.  A  set  {7rn  :  n  £  Nq}  satisfying 
both  of  the  following  conditions  is  the  collection  of  u- 
orthogonal  polynomials: 

1)  Each  7rn  is  a  multivariate  polynomial. 

2)  (7rn,  TTm)^  :=  /  7Tn7Tmda;  =  0  if  n  /  m. 

Jn 

Some  common  collections  of  univariate  orthogonal  poly¬ 
nomials  and  their  associated  weight  functions  arc  given 
in  Table  I  [1], 


Name 

UJ 

n 

Chebyshev  (First  Kind) 

1/Vl  —  X2 

[-Li] 

Chebyshev  (Second  Kind) 

Jl  —  X2 

[-i,i] 

Legendre 

1 

[-i.i] 

Laguerre 

e~x 

[0,  oo ) 

Hermite 

e  1 

R 

TABFE  I:  Some  common  univariate  orthogonal  poly¬ 
nomials,  their  associated  weight  functions  uj,  and  their 
domains  <2. 

Fet  {vrn  :  n  £  Nq}  be  a  collection  of  multivariate 
orthogonal  polynomials  relative  to  the  weight  uj.  Suppose 
the  orthogonal  polynomials  have  been  normalized  so  that 
1 1 7rn  1 1  =  1  for  each  n.  Since  the  orthogonal  polynomials 
relative  to  the  weight  uj  form  a  complete  orthogonal  basis 
of  any  function  /  £  L^(Wid)  can  be  expressed 

as  _ 

/  =  X]  Cn7Tn'  (2) 

where  the  equality  is  understood  to  mean  convergence  in 
norm  and  the  Fourier-like  coefficients  are  given  exactly 
by 

Cn=(/,7Tn)w  (3) 

The  expression  in  Equation  (2)  is  called  the  exact  gen¬ 
eralized  Fourier  series  of  /. 

The  computation  of  Equation  (2)  presents  several  dif¬ 
ficulties.  Not  only  does  the  full-grid  generalized  Fourier 
series  involve  infinitely  many  terms,  but  the  Fourier- 
like  coefficients  can  be  difficult  or  impossible  to  com¬ 
pute  exactly  even  for  seemingly  innocuous  functions 
/.  Moreover,  it  is  challenging  to  implement  numerical 
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integration  techniques  achieving  the  desired  accuracy 
due  to  the  highly  oscillatory  behavior  of  the  integrands. 
Therefore  methods  to  estimate  the  generalized  Fourier 
series  (2)  must  tackle  both  the  series  truncation  and 
coefficient  computation  issues  while  ensuring  that  the 
approximation  error  for  one  does  not  overwhelm  the 
other. 

To  tackle  the  series  truncation  difficulty,  we  explore 
three  sets  of  multi-indices.  Fix  a  positive  integer  N  and 
consider  the  sets  of  candidate  full-grid  multi-indices 

Y$  =  {n  £  Ng  :  n  <  iv}  , 

Tjy  =  |n  G  Nq  :  JlnlU  =  Tv}  , 

and  the  sparse-grid  hyperbolic  cross-shaped  multi¬ 
indices 

&N  =  jn.Nj:n(,  +  l)a  +  l| 

Clearly  \Y$\  =  0(Nd )  and  \S%\  =  0(N  log'7-1  N).  It 
can  be  shown  that  Sff  C  Yd;  [3].  These  pre-determined 
multi-indices  arc  commonly  used  in  spectral  and  pseu- 
dospectral  methods  [9],  [11],  [15],  [16]. 

Define  the  full  and  sparse-grid  truncated  generalized 
Fourier  series  as 

fw  ■=  Y  Cn7Fn  (4) 

new 

for  W€{Y*,T*,S*}. 

The  error  rates  for  the  truncated  orthogonal  polyno¬ 
mial  representations  of  a  function  (4)  depend  on  the  order 
of  the  approximating  polynomials  as  well  as  the  behavior 
of  the  function  itself. 

Theorem  2.1:  For  any  s  >  0  and  /  £  L^(Wl), 

\\f-fw\\u<N-8\\f\\«; 

where  W  £  { Y^,  Tfj.  5'^}  and  ||  •  ||K.  is  the  Korobov 
norm  of  order  s  induced  by  the  inner  product 

d 

(, f,9)Ks  =  Y  (f^njuigiTTnjuYlil+rij)23. 

nSNg  j= 1 

The  proof  of  Theorem  2. 1  regarding  the  accuracy  of 
the  full  and  sparse-grid  truncated  generalized  Fourier 
series  appears  in  [9],  [11]  for  Hermite  polynomial  bases. 
The  proof  can  easily  be  generalized  to  bases  of  other 
orthogonal  polynomial  classes. 


B.  Fourier-like  coefficients 

To  tackle  the  coefficient  computation  difficulty,  an 
optimization  problem  with  a  collocation  method  will 
be  employed.  To  this  end.  Let  A  C  be  a  finite 
collection  of  nodes  with  |A|  =  m,  enumerated  so 
that  A  =  {xi,X2, . . .  ,xm}.  Choose  N,  let  W  £ 
{Yy,  Tfj,  S(ij}  and  say  \W\  =  p.  Enumerate  W  so 
that  1Y  =  {ni,  ri2, . . . ,  np}.  Define  the  jth-row  kth- 
column  entry  of  the  collocation  matrix  X  £  Mmxp  by 
Xjtk  =  7 Tfc(xj),  and  let  the  jth  entry  of  the  vector  f  £  Mm 
be  given  by  f,  = 

The  collocation  method  seeks  a  sparse  collection  of 
coefficients  {cn  :  n  £  W}  that  satisfies 

/ (xj)  =  CnTTn(Xj)  Vxj  £  A.  (5) 

new 

More  succinctly,  Equation  (5)  can  be  expressed  as 
f  =  Ac. 

Although  the  Jacobi-like  collocation  matrix  X  is 
full  rank,  it  is  not  necessarily  square.  Indeed,  we  arc 
interested  in  the  case  m  <C  p.  One  could  relax  the 
equality  in  (5)  and  seek  c  that  minimizes  ||f  —  Xc||2- 
This  is  exactly  a  least  squares  optimization  problem  with 
solution  c  =  where  X t  is  the  pseudoinverse  of  X. 
However,  the  least  squares  solution  does  not  promote 
sparsity.  Unless  one  knew  a  priori  the  exact  sparse  set 
IY  to  use,  it  is  unlikely  using  least  squares  would  result 
in  a  sparse  representation  of  /. 

Instead,  we  propose  approximating  the  collocation  co¬ 
efficients  using  the  Dantzig  selector,  which  is  a  solution 
to  the  optimization  problem 

minimize||c||i 

,  T  II  ,  (6) 

subject  to  || D  lXT {Xc  —  f ) || ^  <  5 

where  D  is  the  rn  x  rn  diagonal  matrix  normalizing  the 
columns  of  X  and  the  scalar  <5  >  0  is  small.  Let  c  be  a 
solution  of  the  optimization  problem  (6).  Then  a  sparse 
generalized  Fourier-series  approximation  to  the  function 
/  is  given  by 

f  =  Y  71-11  •  (?) 

new 

Several  methods  exist  to  quickly  and  accurately  solve 
the  optimization  problem  (6),  including  an  alternating 
direction  method  [10]  and  an  iterative  method  based 
upon  proximity  operators  [12], 

Before  performing  the  error  analysis  of  the  approxi¬ 
mation  (7),  some  notation  must  first  be  presented.  For  a 
domain  fl  C  M'7,  let  C(Q)  denote  the  class  of  all  real¬ 
valued  continuous  functions  on  O,  and  let 

Cm(n)  =  {/  :  /(“)  £  C,  V|a|  <  m}. 
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For  any  positive  integer  m,  define  a  seminorm 

I  '  |xm(n)  on  Cm(Q)  as 

\f\xm(Ci)  =  max  j|/(o)(x)|  :  x  <E  Q,  |a|  <  m|  . 

For  any  function  /  :  -»  R,  let  the  supnorm  of  / 

restricted  to  the  domain  ft  be  given  by  ||/||oo,n- 

The  eiTor  in  the  sparse  orthogonal  polynomial  inter¬ 
polation  of  a  function  /  greatly  depends  on  the  highest 
order  of  the  polynomials  used  in  the  summation  as  well 
as  the  smoothness  and  behaviour  of  /  itself. 

Theorem  2.2:  If  the  nodes  A  = 
{x  :  7TM+i(xj)  =  0}  C  are  a  rectangular  grid 

of  the  zeros  of  the  M  +  1th  orthgonal  univariate 
polynomial  and  if  /  £  L2  (Md)  n  (7m(Md),  then  for 
any  simply  connected  domain  fl  C  Md  there  exists  a 
constant  c  >  0  such  that 


/-  E 

ne« 


Cn7Tr 


<  c2 


—mM 


I/I 


*E(n) 


oo,Q 


The  proof  of  Theorem  2.2  appears  in  [11]  for  multivari¬ 
ate  Hermite  orthogonal  polynomials.  The  generalization 
to  other  classes  of  orthogonal  polynomials  is  straightfor¬ 
ward. 

Finally,  since  the  polynomials  {7rn}  are  orthonormal 
relative  to  w,  the  error  in  approximating  /  by  the  in¬ 
terpolated  polynomial  /  =  ninety  Affin  is  characterized 
by  the  difference  in  their  coefficients: 


2 

/-/ 

UJ 

^  '  (dn  Cn  j  7T n 

new 


=  Cn  -  c 


2 

n|lu>  > 


which  is  guaranteed  to  be  small  provided  c  has  a  small  t\ 
norm.  Indeed,  the  constrant  in  (6)  forces  the  coefficient 
error  to  be  small  between  the  recovered  polynomial  and 
the  interpolated  polynomial. 


III.  Numerical  Examples 

In  the  following  experiments,  the  sparse  generalized 
Fourier  series  approximation  is  computed  for  various  test 
functions  using  the  collocation  approach  presented  in 
Section  2.  In  Experiment  1,  various  node  and  multi-index 
selection  methods  are  used  for  several  test  functions. 
The  accuracy  of  the  approximation  for  each  test  function 
is  presented.  To  illustrate  the  utility  of  the  proposed 
method,  the  collocation  method  presented  above  is  used 
to  quickly  generate  the  Gausian-Hermite  moments  of 
several  images  in  a  training  set  in  Experiment  2,  which 
are  then  used  for  classification  of  rotated  images  in  a 
testing  set. 


The  Hermite  polynomials  arc  used  as  the  orthogonal 
basis  in  the  sparse  representations.  Consider  the  collec¬ 
tion  of  orthogonal  univariate  Hermite  polynomials  with 
associated  weight  function  uj(x.  y)  =  The  first 

few  univariate  Hermite  polynomials  prior  to  normaliza¬ 
tion  are  ttq(x)  =  1,  7Ti(x)  =  2x ,  7r2(x)  =  4x2  —  2,  and 
satisfy  the  three  term  recurrence 

7Tn+l(x)  =  2x'7Tn(x)  -  2n7Tn_l(x) 

for  n  (=  N.  Since  oj  is  separable,  the  multivariate  Hermite 
polynomials  arc  products  of  the  univariate  ones: 

7Tn(x)  =  7Tni(xi)7rn2(x2)  •  ■■TTnd(xd). 

The  polynomials  arc  then  normalized  as 

TTn  F-  7Tn/ \J n!2nV/7T. 

To  find  the  solution  to  the  optimization  problem  (6), 
we  use  the  proximity  operator  based  iterative  method 
recently  proposed  in  [12].  The  proximity  operator  based 
method  is  straightforward  and  uses  only  two  soft  thresh¬ 
olding  operations  at  each  iteration,  requiring  0{Arnp) 
multiplications  in  each  iteration  for  a  m  x  p  collocation 
matrix  X. 

All  experiments  are  performed  in  MATLAB  2014a  on 
a  PC  with  an  Intel  Core  i7-3630QM  2.40  GHz  processor 
and  16GB  RAM  running  Windows  7  Enterprise  with 
machine  precision  2.2204e  —  16. 

Experiment  1. 

In  this  experiment,  the  sparse  generalized  Fourier  se¬ 
ries  approximation  several  given  functions  are  computed. 
Suppose  d  =  2  and  consider  the  functions 

h{x,y)  =  x2y 2 
h{x,y)  =  xV 
f3(x,y)  =  xey/2. 

It  is  easy  to  check  that  fj  <G  L2  (M2)  for  each  j. 

To  create  the  collocation  matrix  X  for  use  in  the 
scheme  (6),  we  use  the  full-grid  rectangular,  triangular 
and  a  sparse-grid  hyperbolic  cross  shaped  set  of  can¬ 
didate  multi-indices,  and  Sfj,  along  with  nodes 

formed  by  a  rectangular  collection  the  zeros  of  the  M  de¬ 
gree  univariate  Hermite  polynomials.  The  parameter  N 
determines  the  highest  degree  hermite  polynomial  used 
in  the  approximation,  and  the  parameter  M  determines 
the  number  of  nodes  used  as  collocation  points. 

Figure  1  displays  the  sparse  generalized  Fourier  Her¬ 
mite  approximations  of  the  functions  /2,/ 3  with  the 
coefficients  computed  as  in  (6),  and  Table  II  gives 
the  error  in  computing  the  coefficients  as  in  (6)  of 
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Fig.  1:  The  exact  function  (left),  the  approximated  generalized  Fourier  series  (center)  and  the  pointwise  approxi¬ 
mation  error  (right)  for  the  functions  fi  (top  row)  and  / 3  (bottom  row)  in  Example  1.  The  approximations  for  fi 
use  Y52  as  the  collection  of  multi-indices  and  M  =  5.  The  approximations  for  fy  use  Y±0  and  M  =  10. 


w 

M 

2 

3 

4 

N 

5 

6 

7 

8 

9 

N-l 

3.1250e-01 

1.8750e-01 

1.5830e+00 

6.9389e-17 

1.3597e-16 

2.3967e-16 

8.7694e-16 

2.0015e-16 

N 

1.8750e-01 

2.3097e-16 

6.9389e-17 

1.6303e-16 

2.3967e-16 

1.2178e-16 

2.0015e-16 

1.2478e-15 

N+l 

2.3097e-16 

6.9389e-17 

1.6303e-16 

2.3967e-16 

1.2178e-16 

2.0015e-16 

1.2478e-15 

5.4918e-16 

N-l 

3.0619e-01 

1.7678e-01 

1.7922e-01 

6.9389e-17 

7.9722e-17 

6.2450e-17 

3.0938e-16 

6.7987e-17 

<e^ 

N 

1.7678e-01 

3.0619e-01 

6.9389e-17 

1.2432e-16 

6.2450e-17 

9.3095e-17 

6.7987e-17 

7.8505e-17 

N-l 

3.0619e-01 

1.0607e+00 

1.2432e-16 

1.6997e-16 

9.3095e-17 

4.6548e-16 

7.8505e-17 

5.1962e-15 

N-l 

3.0619e-01 

1.7678e-01 

3.0619e-01 

1.0607e+00 

2.31 17e+00 

4.0620e+00 

1.1938e-16 

1.6464e-15 

<n£; 

CO 

N 

1.7678e-01 

3.0619e-01 

1.0607e+00 

2.31 17e+00 

4.0620e+00 

6.3122e+00 

1.6464e-15 

2.5542e-15 

N+l 

3.0619e-01 

1.0607e+00 

2.3117e+00 

4.0620e+00 

6.3122e+00 

9.0623e+00 

2.5542e-15 

6.4896e-15 

N-l 

7.2220e-01 

7.3931e-01 

8.2504e-02 

1.3435e-02 

1.6915e-03 

2.6417e-04 

2.0687e-04 

2.065  le-04 

>? 

N 

9.5017e-02 

1.6675e-02 

2.1038e-03 

2.5979e-04 

2.0466e-04 

2.0642e-04 

2.0650e-04 

2.065  le-04 

N+l 

1.0174e-02 

1.2818e-03 

1.9914e-04 

2.042  le-04 

2.0642e-04 

2.0650e-04 

2.065  le-04 

2.0627e-04 

N-l 

7.1773e-01 

6.6775e-01 

1.6675e-02 

2.1038e-03 

2.5979e-04 

2.0466e-04 

2.0642e-04 

2.0650e-04 

N 

5.0888e-02 

1.0174e-02 

1.2818e-03 

1.9914e-04 

2.042  le-04 

2.0642e-04 

2.0650e-04 

2.065  le-04 

N+l 

1.4089e-01 

4.8298e-02 

1.5351e-02 

3.8721e-03 

9.1676e-04 

3.1805e-04 

2.1067e-04 

2.0025e-04 

N-l 

6.4201e-01 

6.6292e-01 

1.4089e-01 

4.8298e-02 

1.0762e-01 

4.4839e-02 

9.2400e-02 

2.5983e-02 

(M  ^ 

Co 

N 

1.1713e-02 

1.4089e-01 

3.2214e-01 

1.0762e-01 

1.7750e-01 

9.2400e-02 

1.6179e-01 

4.6104e-02 

N+l 

1.4085e-01 

3.2214e-01 

5.3508e-01 

1.7750e-01 

2.6136e-01 

1.6179e-01 

2.5687e-01 

6.7726e-02 

TABLE  II:  The  £2  vector  norm  error  1 1 c  —  c || 2  for  the  generalized  sparse  Fourier  series  expansions  f\  and  fy  from 
Experiment  1  for  various  combinations  of  N  and  M. 
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the  test  function  f\  for  each  multi-index  scheme  with 
M  €  {N  -  1,N,N  +  1}  for  TV  =  2, 3,...,  9.  In 
particular',  the  exact  Fourier  Hermite  series  expansion 
of  f\  has  exactly  4  nonzero  coefficients  located  at  the 
multi-indices  I\  =  {(0,  0),  (0,  2),  (2,  0),  (2,  2)}.  As  seen 
in  Table  II,  the  locations  and  values  of  the  nonzero 
coefficients  were  approximated  well  provided  I  \  C  W, 
and  not  well  when  I \  <£_  W.  Similarly,  the  exact  Fourier 
Hermite  series  expansion  of  fa  has  exactly  9  nonzero 
coefficients  located  at  Ij  =  {{ni^nf)  :  n3  =  0,2,4}.  On 
the  other  hand,  the  exact  Fourier  Hermite  series  expan¬ 
sion  of  fy  contains  infinitely  many  nonzero  coefficients. 
However,  as  in  the  Taylor  expansion  of  the  exponential 
function,  the  coefficients  decay  quickly  and  an  accurate 
approximation  can  be  obtain  by  using  only  a  few  terms, 
as  seen  in  Figure  1. 

From  Table  II  one  can  see  the  drawbacks  of  using  the 
predetermined  sparse  multi-index  set  Sfj.  The  parameter 
N  must  be  much  larger  to  achieve  accuracy  on  par 
with  the  full-grid  rectangular  and  triangular  index  set 
implementations,  even  for  a  function  as  simple  as  f\. 

Experiment  2.  In  this  experiment,  the  sparse  Fourier- 
Hermite  series  approximation  is  used  to  perform  pattern 
recognition  of  unknown  images  with  possible  rotation. 
The  Fourier-like  coefficients  are  used  to  compute  the 
rotation  invariant  Gaussian-Hermite  moments  of  orders 
2,3,  and  4,  which  are  then  used  to  classify  the  unknown 
images. 

Consider  a  collection  of  images  in  the  training  set 
Tr  =  {/1, 12, .  ■  . ,  Ik},  where  each  M  x  M  image  has 
been  transformed  to  a  M2  x  1  vector.  For  each  vector  I* 
in  the  training  set,  the  Fourier-like  coefficient  vector  c * 
is  computed  as  in  Equation  (6)  where  the  collocation 
matrix  X  is  slightly  changed.  Instead  of  using  the 
Hermite  polynomials,  each  entry  in  X  is  an  evaluation 
of  a  Hermite  function  Xj^  =  7Tfc(xj)  exp(— 1 1  x^- 1 1 2 /2)  - 
Throughout  this  example,  the  triangular  multi-index  set 
is  used.  Although  the  pixels  in  each  image  naturally 
form  an  equally  spaced  collection  of  nodes,  for  compu¬ 
tational  convenience  the  nodes  are  mapped  to  Zm  x  Zm, 
where  Zm  are  the  zeros  of  the  Mth  Hermite  univariate 
polynomial. 

From  the  Fourier-like  coefficients  approximated  using 
the  optimization  and  collocation  method  as  in  (6),  one 
can  approximate  the  Gaussian-Hermite  moments  of  an 
image. 

Theorem  3.1:  Suppose  /  GE  L2(M2)  and  let  c  be  a 
solution  to  (6)  using  the  orthogonal  hermite  function 
basis.  Then  the  elements  of  c  approximate  the  geometric 


Gaussian-Hermite  moments  defined  by 

f{x,y)'irn{x,y)e~<'x2+y2)/'2  dxdy. 


mn  :  = 


M2 


Proof:  Using  an  orthonormal  hermite  function  basis 
in  (6)  yields  the  approximation 

f{x,y)  =  ^  cnTTn(x,y)e~{x2+y2)/2 
new 


similar  to  (7).  Substituting  the  above  result  into  the 
formula  for  the  geometric  Gaussian-Hermite  moments 
yields 


mn 


f(x,y)Trn(x,y)e  (*2+r)/2  dxdy 
y,  Ct7rt(x,  2/)7T n(x,  y)e~('x2+y~)  dxdy 


'y<  Ct  (vrt,  T^n)u  , 
tew 


which  reduces  to  just  the  element  cn  since  the  hermite 
polynomials  are  orthonormal  relative  to  the  weight  to. 

□ 


Although  Theorem  3.1  is  written  in  the  continuous 
case,  the  discrete  case  follows  by  a  straightforward 
substitution  of  summations  for  integrations.  Using  the 
Hermite  functions  instead  of  the  Hermite  polynomials  is 
important  and  guarantees  that  the  Fourier-like  coefficient 
vector  is  a  close  approximation  to  the  corresponding 
Gaussian-Hermite  moment  for  each  index  in  R2n.  From 
a  small  collection  of  these  approximated  Gaussian- 
Hermite  moments,  one  can  determine  the  rotation  in¬ 
variant  Gaussian-Hermite  moments  fj  of  orders  2,  3, 
and  4.  A  list  of  the  formulas  for  these  rotation  invariant 
moments  is  given  in  the  Appendix. 

For  this  experiment,  seven  images  form  the  training 
set.  These  images  are  8-bit  grayscale  50  x  50  images 
of  simplified  Chinese  characters,  which  were  studied 
in  [17].  The  training  set  is  displayed  in  the  top  row  of 
Figure  2.  Note  that  Images  1-3  are  visually  similar,  as 
are  Images  4-5  and  Images  6-7.  The  testing  set  Ts  is 
formed  by  rotating  the  images  in  the  training  set  by 
0°,  45°,  90°,  135°, . . . ,  315°  followed  by  vectorization. 
For  any  vector  x  in  the  training  or  testing  set,  let  4>x  <G 
M11  be  the  vector  of  the  rotation  invariant  Gaussian- 
Hermite  moments  of  orders  2  thru  4  corresponding  to 
the  vector  x. 

To  perform  classification  of  the  rotated  images,  say 
x  G  Ts  is  classified  as  a  rotation  of  image  y  G  Tr  iff 

||$x  -  ^ylli  <  ||$x  -  ‘Mi  Vz  G  Tr. 
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This  method  of  using  the  rotation  invariant  Gaussian- 
Hermite  moments  of  orders  2,3  and  4  as  the  feature 
vector  yielded  100%  correct  classification  rate  of  all  56 
images  in  the  testing  set.  Values  of  |4>x  —  4>y||i  taking 
x  over  all  7  images  rotated  90°  and  y  over  all  elements 
in  the  training  set  are  shown  in  Table  III.  It  is  clear 
from  the  table  that  although  many  of  the  images  are 
visually  similar,  the  rotation  invariant  Gaussian-Hermite 
moments  are  similar  only  for  rotations  of  the  same 
images.  Of  course,  if  the  Gaussian-Hermite  moments  are 
approximated  well,  then  this  method  is  expected  to  yield 
100%  classification  accuracy.  The  novelty  comes  from 
computing  the  method  used  to  compute  the  Gaussian- 
Hermite  moments.  Employing  the  solution  of  (6)  is  an 
elegant  way  to  compute  the  moments  without  using 
quadrature  schemes  or  hierarchical  methods  [8],  [14], 

IV.  Summary 

In  this  paper,  a  method  for  computing  a  sparse  gener¬ 
alized  Fourier  series  with  orthogonal  polynomial  bases 
of  multivariate  functions  was  presented.  The  method  first 
uses  a  truncated  set  of  multi-indices,  then  approximates 
the  coefficients  using  the  solution  of  an  optimization 
problem  involving  a  collocation  model.  Several  examples 
were  presented  to  illustrate  the  accuracy  and  utility  of 
the  proposed  method.  In  the  first  set  of  experiments,  the 
sparse  generalized  Fourier  series  approximation  of  sev¬ 
eral  test  functions  were  computed,  and  the  approximation 
errors  were  presented.  In  the  second  set  of  experiments, 
the  coefficients  in  the  sparse  generalized  Fourier  series 
were  computed  and  used  as  rotation  invariant  feature 
vectors  to  perform  pattern  recognition  on  unknown  ro¬ 
tated  images. 
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Fig.  2:  The  images  forming  the  training  set  and  testing  set  from  Example  2.  The  top  row  contains  the  seven  images  in 
the  training  set  and  the  remaining  images  form  the  testing  set  of  images  rotated  counterclockwise  0°,  45°,  90°,  135°, 
etc. 


k\j 

1 

2 

3 

4 

5 

6 

7 

1 

6.9480e-09 

4.9948 

4.5607 

7.4822 

10.5865 

6.2996 

8.6762 

2 

4.9948 

2.1222e-09 

2.9397 

8.0141 

11.3291 

7.4386 

7.0308 

3 

4.5607 

2.9397 

3.4982e-09 

6.7547 

10.4790 

6.5229 

7.5554 

4 

7.4822 

8.0141 

6.7547 

9.8842e-10 

6.7450 

4.0142 

5.3556 

5 

10.5866 

11.3291 

10.4790 

6.7450 

9.2692e-10 

6.1498 

6.7113 

6 

6.2996 

7.4386 

6.5229 

4.0142 

6.1498 

3.9752e-09 

3.5455 

7 

8.6762 

7.0308 

7.5554 

5.3556 

6.7113 

3.5455 

1.9313e-09 

TABLE  III:  Values  of  ||$T  —  $/fc||l5  where  x  is  the  element  7y  <E  Tr  with  a  90°  rotation. 
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Appendix 


Below  are  the  complete  and  independent  set  of  rota¬ 
tion  invariants  of  Gaussian-Hermite  moments  of  orders 
2,  3  and  4  [18]. 

0i  =777-20  +  m-02 

fa  =(m30  +  mu)2  +  (m03  +  m2i)2 

fa  ={rn2 o  -  mo2)[(m3o  +  mi2)2  -  ( m03  +  m2i)2} 

+  4mii(m30  +  mi2)(m03  +  m2 1) 

fa  =mn[(m30  +  mfa2  -  (m0 3  +  77i2i)2] 

-  (m20  -  mo2){m3o  +  mi2)(m03  +  m2i) 
fa  =(,m30  -  3mi2)(m30  +  mi2)x 

x  [fan 30  +  mi2)2  -  3(m03  +  m2i)2] 

+  (m03  -  3m2i)(m03  +  m2i)x 
x  [(m0 3  +  m2i)2  -  3(m30  +  mi2)2] 
fa  ={m30  -  3mi2)(m03  +  m2i)x 

x  [(m03  +  m2i)2  -  3 (m30  +  mi2)2] 

-  (3m2i  -  nio3)(m3o  +  mi2)x 

x  [(m30  +  mi2)2  -  3(m03  +  m2J 

fa  =77740  +  2m22  +  777-04 

fa  =(w-40  -  ^04)[(w-30  +  TO12)2  -  (m2i  +  77703)2] 

+  4(77731  +  777-13)  (77730  +  777i2)(?772i  +  777-03) 
fa  =(?773i  +  777i3)  [(777-30  +  777i2)2  -  (?772i  +  77703)2] 

010  =( 777-40  -  6?7722  +  77704)  [(77730  +  777l2)4 

-  6(77730  +  777-l2)2(?772l  +  771-03  +  (^7721  +  77703)4] 

+  16(777-31  -  777i3)  (77730  +  777,12)  (777-21  +  77703)  X 
X  [(77730  +  777-l2 )2  -  (77721  +  77703)2] 

011  =(777-40  -  6?7722  +  77704)  (777.30  +  777l2)(?772l  +  ?7703)X 
X  [(m2i  +  777.03)2  -  (77730  +  777l2)2] 

-  (77731  -  777i3)  [(77730  +  7?7i2)4 

-  6(77730  +  777l2)2(r?72l  +  77703)2  +  (77721  +  77703)4]- 


